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Abstract 

We analyze the power of quantum adiabatic evolution algorithms 
(QAEA) for solving random NP-hard optimization problems within a 
theoretical framework based on the random matrix theory (RMT). We 
present two types of the driven RMT models. In the first model, the 
driving Hamiltonian is represented by Brownian motion in the matrix 
space. We use the Brownian motion model to obtain a description 
of multiple avoided crossing phenomena. We show that the failure 
mechanism of the QAA is due to the interaction of the ground state 
with the ” cloud” formed by all the excited states, confirming that in 
the driven RMT models, the Landau- Zener mechanism of dissipation 
is not important. We show that the QAEA has a finite probability 
of success in a certain range of parameters, implying the polynomial 
complexity of the algorithm. The second model corresponds to the 
standard QAEA with the problem Hamiltonian taken from the Gaus- 
sian Unitary RMT ensemble (GUE). We show that the level dynamics 
in this model can be mapped onto the dynamics in the Brownian mo- 
tion model. However, the driven RMT model always leads to the 
exponential complexity of the algorithm due to the presence of the 
long-range intertemporal correlations of the eigenvalues. Our results 
indicate that the weakness of effective transitions is the leading effect 
that can make the Markovian type QAEA successful. 
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Introduction 


The quantum adiabatic algorithms (QAA) are designed for solving combina- 
torial search and optimization problems based on the quantum evolution of 
the relevant systems [1]. These algorithms are based on adiabatic theory of 
Quantum Mechanics, stating that the quantum state is closely following an 
instantaneous ground state of a slow- varying in time control Hamiltonian, 
provided that the initial state was chosen to be a ground state of the initial 
Hamiltonian. At the initial moment of time the control Hamiltonian has 
a simple form with the known ground state that is easy to prepare and at 
the final moment of time it coincides with the “problem” Hamiltonian Hp 
which ground state encodes the solution of classical optimization problem in 
question. It also can be chosen to reflect the bit-structure and cost spectrum 
of the problem. For example, 

H p = Y j E z |z)(z| (1) 
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Here E z is a cost function defined on a set of 2 n binary strings z = {z Xj . . . , z n } 
Z j = 0, 1, each containing n bits. The summation in (1) is over 2 n states |z) 
forming the computational basis of a quantum computer with n qubits. State 
| Zj)j of the j-th qubit is an eigenstate of the Pauli matrix b z with eigenvalue 
1 — 2 Zj ±1). If at the end of QAA the quantum state is sufficiently close to 
the ground state of Hp then the solution to the optimization problem can be 
retrieved by the measurement. 

It was shown [2, 32] that the system can be trapped during the QAA in 
the local minimum of a cost function for a time that grows exponentially in 
the problem size n. It was also shown [2] that an exponential delay time in 
quantum adiabatic algorithm can be interpreted as a quantum-mechanical 
tunneling of an auxiliary large spin system. 

The paper [4] demonstrates that tunneling can be avoided by a suit- 
able modification of the QAA evolution ” paths” with the same initial and 
final wave functions. The authors presented a general approach using the 
” stochastic” paths and numerical results indicated that the HWP may be 
solved in polynomial time with finite probability. This approach was fur- 
ther developed in [5], where the complete analytic characterization of the 
’’stochastic” path evolution and probability of success was obtained. 
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In present paper, we analyze the performance of the Quantum Adiabatic 
Algorithm (QAA) in the framework of the Random Matrix Theory (RMT). 
We show that the runtime of the algorithm should be typically exponentially 
large in the number of qubits n.in order to have a polynomial probability 
of success for the quantum evolution algorithm. This implies that the QAA 
typically has exponential complexity. 

2 Quantum Adiabatic Evolution Algorithm 

In a standard QAA [1], one specifies the time-dependent Hamiltonian H(t) 

H(t) = a(t) H d + (3 (t) H P , (2) 

with the initial and terminal conditions 

a(0) = 1, 8 (0) = 0, (3) 

a(T) = 0, P(T) = 1. 

The Hamiltonian (2) guides the quantum evolution of the state vector 
from t = 0 to t = T. the run time of the algorithm, and Hp is the “prob- 
lem” Hamiltonian given in (1). In case when the parametric evolution of 
the Hamiltonian (2) is sufficiently slow, the Adiabatic theorem assures that 
the system initially occupying the ground state of the driver Hp = H (0), 
parametrically evolves into the ground state of Hp = H(T). The standard 
estimate of the run time that provides such adiabatic evolution, is well known 
and given by [1] 


T> 1/A 2 , (4) 

where A is the lowest value of the energy gap between the ground and the 
first excited adiabatic states. 

The ” driver’ Hamiltonian Hp is designed to cause the transitions between 
the eigenstates of Hp . According to the above discussion, no prepares the 
initial state of the system j^(0)) to be a ground state of H{ 0) = Hq. It is 
typically constructed assuming no knowledge of the solution of the classical 
optimization problem and related ground state of Hp. In the simplest case 

n-l _ 

H D = -CY J °l IV’(O)) = 2~ n ' 2 jz) , (5) 

j=0 z 
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where a 3 x is a Pauli matrix for j - th qubit and C > 0 is some scaling constant. 
Consider instantaneous eigenstates |^(r)) of H(r) with energies E k (r) ar- 
ranged in non-decreasing order at any value of r € (0, 1) 

H\<pk) = Ek\4>k), k = (6) 

Expanding the non-stationary wave function of the system in terms of the 
instantaneous adiabatic basis (which is complete for any fixed r G (0, 1)), we 
obtain ([6]): 


N 

MO) = J2 Ck (*) exp 

k= 0 

with N = 2 n - 1 complex coefficients { C k (t)} satisfying the normalization 
condition J^, k =o \C k (t ) | 2 = 1 for t G (0,T), and we use the units with h = 1 . 
Note that each quantity \C k (t)| 2 (k = 0, ... N ) is the probability to find the 
system in the instantaneous adiabatic eigenstate \4> k ( t )) at time t. According 
to the above discussion and taking into account the normalization condition, 
the initial conditions for (7) are given by 


-i / dt'E k (t') 


I <t>k (t)), 


( 7 ) 


Co (0) = 1, Cjt (0) = 0, k = l,...,N. (8) 

Under the assumption that the populations of the excited states of the in- 
stantaneous Hamiltonian are small, the Schroedinger equation applied to (7) 
yields [6] 

^ C k (t) = (d> k (t) | ^ |o 0 (t)) exp | -i J dif [E k (?) - E 0 (?)} | , (9) 

describing the dynamics of the non-adiabatic transitions from the ground 
state k = 0 to the states with k = 1 ,N. Taking into account (9), we 
obtain the probability Po (?) for the system to remain in adiabatic ground 
state in the form 


P 0 (t) = l-T (t ) , (10) 

with P (?) being a total probability of escape from the ground state, or failure 
of the QAA, defined by 
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In (11), (Jj) = (d>k (t) \-§i\4>o (t)) is the matrix element of the operator Jj 
and u>k o (t) = Ek (£) — Eo (<) is. the transition frequency between the ground 
and the fc-th excited adiabatic state. In what follows, we will analyze the 
probability (11) for a particular case of a driven RMT model for the problem 
Hamiltonian. 


3 Random Matrix Analysis of QAA 

According to the above discussion, the problem Hamiltonian of QAA should 
have a complex ground state encoding the optimal solution of the initial 
optimization problem. As a simple ” limiting case” of such a complex Hamil 1 
tonian, one may suggest a random Hamiltonian sampled from one of the 
symmetry ensembles of the Random Matrix Theory (RMT) [7]. It has been 
shown that the "random” Hamiltonians belonging to a certain symmetry 
class (say Unitary Ensemble, or GUE), give a good description of the com- 
plex physical systems, which do not possess any specific conservation laws 
and are only characterized by the symmetry class with respect to the change 
of basis (unitary transformations) [7], [8]. We will give more motivations for 
such choice in a subsequent section. In this section, we will analyze a specific 
case of such a "random” problem Hamiltonian,. when the matrix elements of 
the problem Hamiltonian evolve according to the Ornstein-Uhlenbeck pro- 
cesses [18] 


d JM± = - 1 H kj (t) + a^ j , ( 12 ) 

where 7 > 0 describes the dissipation, a > 0 the noise amplitude and z kj ( t ) 
correspond to N 2 independent Langevin sources 

<&; (t ) , br ({)) = Skk'SjfS (t - t ') . (13) 

In (13), Ski stands for the Kronecker symbol and 6 (t — t r ) is a Dirac’s delta 
function. The formal solution of the stochastic differential equation (12) is 
given by 
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( 14 ) 


H kj ( t ) = H kj (0) exp [-7 1 ] + cr [ dt'£ kj (t 1 ) exp [-7 ( t - £')] , 

JO 

where H k j (0) is an arbitrary Hermitian matrix representing the ” initial con- 
dition” for the driven random matrix H k j (t = 0) = H kj (0). 

3.1 Angular Decomposition in GUE 

The ’’angular coordinates” for the Hermitian matrices H belonging to the 
General Unitary Ensemble (GUE) [7] are defined by 

H=UAU+, (15) 

where the diagonal matrix Ay — Sij\ represents the eigenvalues {A^} of H 
which are invariant with respect to the unitary transformations = ijj^ 
characterizing the j-th component of the i-th eigenvector belonging to the 
eigenvalue number i. Therefore, the unitary matrix U represents a particular 
choice of the basis and the vector A is basis-invariant. Note that the unitarity 
conditions U + U = UU+ — I where I is a unit matrix, correspond to the 
orthonormality and completeness of the basis Differentiating (15) and 
taking into account the unitarity conditions, we obtain 

U+dHU = dA + [dQ, dA] , (16) 

where the generalized ’’angular” motion is defined by 

dQ = U+dU, (17) 

and [...] is a commutator, he. [dQ, dA] = dQdA — dAdQ. Making use of (16), 
we obtain the equations of motion for the eigenvalues and angles Q in the 
’’local” basis with U = /, in the form 


dt 



a 


\ ~ \ 


dQ k j g £ 

” A j - Xk j 


(18) 
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From (15), it follows that 


dftfeO 

dt 





A* (t) - A 0 (t) 


E-k ( t ) — Eo ( t ) = U>kO ( t ) , 


( 19 ) 


establishing the connection between the probability (11) and the angular 
decomposition variables. 


3.2 Optimal Trajectories for Driven RMT 

The first of equation (18) describes the many-particle system in ID un- 
dergoing the Brownian motion and interacting with the long-range pairwise 
potential ~ 1/r [7], [8]. Importantly, the dynamics of the eigenvalues A (t) 
does not depend on the eigenvector dynamics ft (£). Since the velocities 
of the Brownian particles contain the noise terms, each realization of the 
general many-particle trajectory in the phase space is drawn from the ran- 
dom distribution characterized by a particular choice of the Langevin source 
trajectories (t). Therefore, the probability of success of the QAA (11) be- 
comes a functional defined on each of the Langevin source trajectories. Since 
these trajectories are random, this probability Pa is characterized by the dis- 
tribution in the general phase space of the problem, with the ’’distribution 
density functional” defined on the realizations of the Langevin sources and 
therefore on the phase coordinates of the system, in the form [23] 


W [A (t ) , Q (t) } = exp {-S [A (t) , ft (t)]} , (20) 

with the action functional S defined through the Lagrangian function 


S = 



A (if) , A (t') , ft (t') , ft (t') 


( 21 ) 


where A = dA/dt and ft = dQ/dt, and L is the Lagrangian. Following the 
standard procedure, the Lagrangian for the system (18) is given by [22] 
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L = 


2a 2 


E 

. k . 


Afc - Ufc ) + cr 2 ^ 


+ ^ (A*; — Aj)“ ( Qaj i [ , 
kj 


( 22 ) 




= M 


J^k 


G 

A fc ~ A , / ’ 


with t'( (A) = dvk (A) /3Afc. Making use of (11), (18), and (19), we obtain 
the expected probability of failure (F (£)) = 1 — (Pq (£)) in the form 


<r (t)) = 


J f DADQW [A, 12] EL /o df J[ dt" 12 fc0 (f) 12 0fc (£") exp {-z$ fc0 } 


/ / DADQW [A, 12] 


(23) 


where (£', £") = / t f <fnj fc0 (r), f DA f DQW [A, 12] denote the functional 
integration over all trajectories with the weight defined as a functional prob- 
ability density (20), and the denominator is included into (23) for normal- 
ization. From (22), it follows that the ’’angular” part enters as an additive 
quadratic term into the Lagrangian, and therefore the integration over the 
angular variables f DQ in (23) can be done explicitly [23], with the result 


f DOW [A, 12] Qko ( t ) Qpk (Q _ 


1 1 DADQW [A, 12] V [A*, (f) — A 0 (£)f 

Combining (23) and (24), we obtain 


5 {t - £') 


(24) 


N 


<r («)> = £ / da 


exp (-S eff [A]) 


dt'- 


fc= i 


/ DA ' exp (—S e ff [A']) J 0 [A* (£') - A 0 (£')] 


r, (25) 


where the effective action S e ff [A] is taken with the ’’truncated” Lagrangian 
Leff 


L *ff = 2^E 


k L 


A k ~ (A) ) + G~v f k (A) 


(26) 



The evaluation of (25) is still very difficult, since it involves the averaging 
over the unknown level distribution. 

For the many-body system with exponential number of particles N — 2", 
the distribution of the system’s trajectories in the phase space is expected to 
be ’’sharp” around the optimal trajectory, defined from the variation principle 
[21]. The variational principle takes the form of minimum action 

S„ in = min {S [A (t), (!(*)]}, (27) 

Given (22). the optimal trajectories are obtained from the minimum action 
principle (27) as the solutions of the Euler-Lagrange equations 


d_ ( dL\ _ fdL\ 
dt \dqJ V%/ 


(28) 


where g* = {A,} represent the generalized coordinates [9j. Combining (22) 
and (28), we obtain 


= a), 


with the effective potential given by 

i 




Combining (22), (29) and (30) yields 


(29) 


(30) 


A ' k (t) = ~U eff ( A), 

OAk 

u«f( A) - 

k 


(31) 


The details of the calculations leading to (31) are presented in Appendix A. 
The effective potential in (31) can also be derived from the Fokker-Plank 
equation describing the time evolution of the N - particle level distribution 
function [18]. In general, that procedure leads to the dynamic equations 
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describing the many-particle system in ID with pairwise interactions. In a 
particular case of the unitary ensemble, the pairwise interaction vanishes and 
the problem reduces to the IV-particle non-interacting system. Therefore, we 
finally have the following description of the optimal trajectories 


A* (t) - 7 2 A* = 0. (32) 

The equations (32) are of the second order and therefore require two initial 
(or boundary) conditions for each k — 0, ...N. Since (18) are the first order 
stochastic differential equations, we only have one initial condition for the 
optimal trajectories 


A* (t = 0) = X k (0) , , (33) 

for each k and the final points A*, (t = T) are not fixed. We impose the 
additional boundary conditions at t — T in the form 


-lX k {T) + 



X k (T) - A; (T) 


= 0 , 


(34) 


which insure that at t = T the system belongs to the equilibrium unitary 
ensemble. Specifically, (34) corresponds to the standard GUE distribution 
[7], [8]. As we show in Appendix B , the boundary conditions (34) can also be 
obtained in the long-time limit for the optimal trajectories in case when the 
terminal conditions at t = T are not fixed and have to be self-consistently 
derived from the optimality property of the trajectories. In fact, we show 
that the terminal conditions for the ” arbitrary 7 ’ optimal trajectories are ex- 
ponentially close to the equilibrium described by (34). This implies that the 
77 spread 57 of the optimal trajectories is exponentially small in the large-time 
limit, confirming the validity of the approach described above. 

The equations (32) with (33) and (34), provide a close set characterizing 
the optimal trajectories of the initial many-body system. In the subsequent 
section, the probability of success of the QAA will be evaluated using the 
solutions for the optimal trajectories. 


3.3 Probability of Success for the QAA 

As we have mentioned above, for the system with exponential number of 
particles A, the distribution of trajectories in the phase space is expected 
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to be ’’sharp” around the optimal trajectory. Following this argument, we 
observe that with exponential accuracy [21] 


<T(t)> 


' V /* n 2 

= Y' / dt' ^ 

[A; (?) - Aq (t')f 


(35) 


where {Aj“} correspond to the optimal trajectories. In order to simplify the 
notations, w'e will denote {At} = {A*,}, i.e. drop the superscript. This will 
not cause the confusion, since we will only consider the optimal trajectories 
in the rest of the paper. From (32) with (33) and (34), we obtain 


A* (f) = A, (0) 


sinh [7 (T — £)] i sinh (yt) 

sinh (7 T) ' k ' sinh (7 T) 


(36) 


The final ’’coordinates” {A^ (T)} are defined from (34) and correspond to the 
standard equihbrium GUE distribution [7], [8] at t = T. Substituting (36) 
into (35), we obtain 


(r (D) 


N T 

= V f 


dt'- 


cr sinh 2 (7 T) 


Jo [tUfco (0) Sinh [7 (T - f)] + u k 0 (T) sinh (7^)] 
The integration in (37) can be done explicitly, with the result 


2' 


(37) 


(F(T)) 
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(7 

7 


N 

sin.h{-fT)^2 
k= 1 


1 

w*:0 (0) cufco (T) 


(38) 


Making use of the terminal condition (34), the probability (38) can be esti- 
mated as 


<r (T)> = sinh ( 7 r) Ao (T) • (39) 

Since the Hamiltonian at t = 0 is given by (5), the average inverse frequency 
in (39) is bounded 


n \&ko (0) 


On the other hand, according to (34), we have 


( 40 ) 
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(41) 


A 0 (T) 


aV2N 

~7T' 


Since the spectrum of the problem Hamiltonian is bounded in the largen 
limit, the ground state energy Ao (T) should not depend on n. This is satisfied 
if 



N’ 


(42) 


which is a common scaling in RMT models [7], As we will see in the subse- 
quent section, the scaling (42) naturally occurs in the QAEA framework. In 
Appendix C, we show that the equilibrium iV-level distribution only depends 
on the parameter a\ = a 2 / 7, confirming that the condition (42) is sufficient 
for the spectrum to be bounded in a finite interval. 

Given the estimate (40), the summation over the energy spectrum is 
performed using the exact equilibrium condition (34), and therefore the result 
(39) is exact. This means that (39) is based on the exact summation over 
the discrete spectrum valid for any total number of states N and does not 
refer to any continuous level distribution arising in the large N limit. In the 
large N limit, the equilibrium density of states (DOS) for the GUE ensemble 
is given by the famous ’’semicircle” law [7], [8] 


2 N / 

P( A ) = ^V A o-A 2 , — A o < A < A 0 , (43) 

with the normalization condition J dXp (A) = N. One should note that 
near vicinity of the ground state, the level distribution functions may have 
large corrections, which may significantly change DOS near the ground state 
in comparison to the bulk distribution [14]. However, one can show that 
the main contribution to the probability (38) comes from the bulk of the 
distribution, where there is a large number of states and (43) is valid in the 
large N limit. Combining (38) and (43), we obtain analogous to (39) 


2A 0 

<r (t» = d sinh S ^ (2A ° - “)■ (44) 

0 

From (44), it follows that the integral over the spectrum in the r.h.s. of (44) is 
accumulated at the upper bound of the integration, since the integrand only 
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has integrable singularity at u> = 0. This indicates that the main contribution 
to (44) and therefore to (38), comes from the upper part of the spectrum. 
After evaluation of the integral, (44) yields 

(r(T) , =smhhT) (^(0 ))y£' (45> 

which reduces to (39) after taking (41) into account. Based on the above 
arguments, we can employ the semicircle distribution to evaluate the average 
spacing This is done in Appendix D, and the probability (39) takes 

the form 


(r(T)H — A o( r)sinh( 7 r). (46) 

7 rn 

Since yT ~ 1 and the absolute value of the ground state energy Ao (T) does 
not depend on n, it follows from (46) that the failure probability of the QAA 
is (T (T)} ~ 1/n, implying that the algorithm can be successful in the largen 
limit. As we will show in the subsequent section, this is a consequence of the 
Markovian nature of the Brownian motion model considered above. In the 
next section, we will discuss this in more details and compare to the results 
in the non-Markovian case. 

The result (46) could not be obtained using a standard estimate (4) with 
A = o-’io being the smallest gap between the ground state and the first 
excited level of the total Hamiltonian (?) As it follows from comparison of 
(38) and (34), the main contribution to the failure probability comes from' 
the bulk of the spectrum, implying that the failure of the QAA occurs due 
to the interaction of the ground state with the ” cloud” formed by all the 
excited states with k = 1, ...N. In other words, it is the large number of 
possible transitions from the ground state that may cause the failure of the 
QAA, not the structure of the minimum gap in the energy spectrum. The 
main reason for this is that the phase factors present in (11) vanish after the 
integration over the angular variables (24). Note that these phase factors 
contain the transition frequencies and lead to the low transition probabilities 
in adiabatic regime. This is not the case in the Brownian motion model 
considered above, when all the excited levels work together and the excitation 
can be effectively viewed as a transition to the ” cloud” of the excited states. 
This is in qualitative agreement with [13], confirming the conclusion that in 
the driven RMT models, the Landau-Zener mechanism of dissipation [12] is 
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not important. Indeed, the theory of dissipation based on the Landau-Zener 
mechanism assumes [ 11 ] that the transitions happen due to the pairwise 
"level crossing'’ events, rather than the ” collective" interactions picture that 
follows from our analysis. 

In the regime when the QAA is successful considered above, the average 
minimal gap is A ~ 1/N, and the formal application of the standard esti- 
mate gives T > N 2 , which is has nothing to do with the result (46). One 
should note that in some cases, the bound (4) can be improved by optimally 
adjusting the control parameter corresponding to the rate of the quantum 
adiabatic evolution [ 1 ]. In our case, this would correspond to introducing 
the time dependent rate parameter 7 (t) in (32), and then optimizing with 
respect to the function 7 (t). However, as we show in Appendix E, the result 
(38) is essentially invariant with respect to the choice of 7 (£), and therefore 
(46) does not depend on the functional form of the rate parameter. The 
main reason is that, according to the above discussion, the QAA fails due 
to the collective interaction of the ground adiabatic state with the cloud of 
excited states, and therefore there no well defined "dangerous" moments of 
time corresponding to the avoiding crossing [ 12 ], [ 11 ]. 

One may attempt to correct the standard estimate (4) using the average 
level spacing A = (u/ 10 ) in the vicinity of the ground state instead of the 
average spacing. The point is that the spacing near vicinity of the edges of 
the spectrum may be significantly enhanced in comparison to the average one 
[14]. According to the above discussion, these estimates can not be correct, 
since they assume that the main contribution to the failure probability comes 
from the low-energy part of the spectrum. 

4 Motivations for Brownian Motion Model 

In this section, we will show that the specification considered in a previous 
section is actually quite general, since it can be mapped onto ( 2 ), with the ap- 
propriate choice of the ’’trajectory” parametrized by the functions a(t) and 
(3 ( t ), and the problem Hamiltonian Hp described by a random matrix belong- 
ing to the Gaussian Unitary ensemble (GUE). In general, such mapping yields 
the time-dependent coefficients { 7 , a} characterizing the Ornstein-Uhlenbeck 
process. We will specify a particular parametrization that leads the constant 
coefficients, and therefore the stationary distributions. 

Consider the QAEA ( 2 ) with a random GUE problem Hamiltonian Hp = 
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V = H ( T ) and the driver Hp = H (0). Since V belongs to the GUE, its 
matrix elements are independent and Normally distributed, V kj N (0, a 0 ). 
Taking this and (2) into account, we derive the probability distribution for 
the matrix elements of the Hamiltonian H(t) in the form 


P [H, t ] = C N wn exp 
kj 


H kj -a(t)H( Q) kj 


2 


(47) 


where Cy (£) is the normalization constant. Suppose there is a change of 
parametrization r = r (t) in (2) that maps the total Hamiltonian dynamics 
(2) onto the Ornstein-Uhlenbeck process (12). The process (12) leads to the 
following probability distribution 


with 


P [H, r j 



[H kj -H kj (r )] 2 ' | 

4 D (r) j ’ 


(48) 


H kj (r) = H (0) kj exp (— yr) , 

D (r) = 1 [ dr a 2 (r) exp [~2j (r - r')] , 


(49) 


and we assumed that the effective variance a 2 (r) may be time-dependent. 
Comparing (47) and (48), we obtain the conditions on the parameters 7 and 
cr (r), in the form 


exp(— yr) = a(t), 



dr cr (r') exp (27/) 


VqP 2 ( t ) exp (2 7t) , 


(50) 


and therefore 
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r(t) = -1 In [<*(*)]> 

a 2 (r) = crlexp(-2yT)^-{0 2 [t(T)]exp(2-/T)} , 


(51) 


implying that the effective variance a 2 (r) is indeed in general time-dependent. 
In what follows, we consider the QAEA (2) 


H(t) = H (0) exp (-7 1 ) + V [1 - exp ( — yt)] , (52) 

where 7 is a real parameter characterizing the speed of the parametric evo- 
lution. The parametrization (52) is a particular case of (2) with 

a(t) = exp (-7 1 ) , 0(t) = 1 - exp (-7 1 ) , T ~ (53) 

7 

and the parameters a(t) and (3 (f) being the monotonic functions for t 6 
[0, Too] of (52). From (51), we derive the parameters of the corresponding 
Ornstein-Uhlenbeck process in the form 


T{t) = t, 

a 2 (t) = 270-“ [1 -exp(-Tf)] , 


(54) 


implying that the trajectories (52) are mapped onto the Ornstein-Uhlenbeck 
process with time-dependent variance parameter. From (54), it follows that 
the time dependence of the effective variance is slow and achieves the "sat- 
uration” during the runtime of the algorithm, when t ~ 1/7. As we show in 
Appendix F, there is a particular class of trajectories corresponding to the 
Ornstein-Uhlenbeck process with constant parameters [17]. 

As it was pointed out in [17], the mapping (51) is exact in a sense that the 
level distribution obtained from the corresponding Ornstein-Uhlenbeck model 
is identical to the level distribution of the parametrically driven Hamiltonian 
(2), at each instant of time. Hov/ever, the actual level dynamics of the two 
models are quite different. For this reason, the level correlations calculated 
at different time moments are not equivalent for (12) and the corresponding 
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(52). In particular, the driven system (52) contains the ’’quenched” disor- 
der terms leading to the intertemporal correlations, which are absent in a 
Markovian system (12), [17]. In case of the QAEA with the random Prob- 
lem Hamiltonian (52), the intertemporal correlators can be easily estimated. 
Making use of this and the knowledge of the effective level dynamics obtained 
by mapping onto (12), we will be able to estimate the probability of success 
for the model (52). 

The probability of failure for the QAA is obtained from (11) in the form 


* rt r t H{ t') H(t ") f ) 

r ' (t) = £ l * l -JoL(r) exp { l (T) } • 

(55) 


where we took into account the identity 


= »(<) 


fcO ^*0 (t) 

with H ( t ) = 4M-. In case (52), we have ^ = —7 ( H — V ), and (55) yields 


(56) 


r(«) = 7 V?E[c 4 (()| ! , 


k = 1 


exp <-if* druj k0 (r)} 

c ‘ (i) - L* ■«.(*) — - 


(57) 


where we took into account the correlations (V k0l V 0k ) = crj and that the 
Hamiltonian H (t) is diagonal in the adiabatic basis. 

The second of (57) represents the transition amplitude, which typically 
contains two components. In the adiabatic limit, the leading contribution to 
the integral (57) comes from the endpoints of the integration and corresponds 
to the corrections to the adiabatic eigenstates due to the finite rate of the 
parametric transformation (52). The second contribution comes from the 
actual transitions between the '’true” dynamic states with the non- adiabatic 
corrections to the states taken into account. The ’’transition” term is ex- 
ponentially small in the adiabatic limit and exhibits the non-analytical de- 
pendence on the adiabaticity parameter ([6]). On the other hand, the ’’state 



correction” term typically has a power-law dependence on the adiabaticity 
parameter. 

The general setting of the QAEA requires that the solution of the initial 
optimization problem is encoded in the adiabatic ground state of the Hamil- 
tonian [1], Therefore, in order to evaluate the success of the QAEA, one 
has to take into account both ’’correction” and ’’transition” contributions to 
the integral (57). The correction term gives the leading contribution and 
therefore determines the success of the QAEA. This distinguishes between 
the approach adopted in the present paper and the one from ([10]). In ([10]), 
the authors evaluated the failure of QAA based on the results ([11]) for the 
dissipation rate in the driven RMT ensembles. The dissipation was estimated 
in ([11]) based on the Landau-Zener theory, which takes only transition com- 
ponents in (57) into account. While this was a reasonable approximation to 
the estimation of the total dissipation rate, this approach is not applicable 
to the analysis of the QAA. The point is that the evaluation of the QAA 
performance requires taking into account any deviations from the adiabatic 
ground state, regardless whether they occur due to the ’’true” transitions to 
the ’’corrected” states, or the corrections to the states themselves. There- 
fore, it is essential that both contributions to (57) be taken into account. 
Moreover, from the above discussion, it follows that the ’’correction” terms 
give the leading contribution, and therefore the ’’transition” effects can be 
neglected in the adiabatic limit. This is consistent with the results of ([13]). 
The ’’correction” contribution in (57) is given by ([6]) 


c k (t) = 


— % 


jl-exp -i J dt'io k0 (t') j. (58) 


[ktfcO ( t )] 2 

Combining (57) and (58), we obtain the failure probability of QAA 


N 


r(T) = 


No PT 


sin" 


dt'uko (0 


(59) 


In order to evaluate the failure probability (59), we have to specify the 
level dynamics of the driven model (52). Making use of the mapping (54) en- 
ables us to apply the optimal trajectory methods developed for the Brownian 
motion model. The optimal trajectories characterized by (36), and we obtain 
for the transition frequencies along the optimal trajectories in the form 
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(60) 


^fco (t) = ^'fco (0) 


sinh [7 (7 1 — 4)] 
sinh (7T) 


+ ^’fco (T) 


sinh (7 1) 
sinh (7 T) ' 


Substituting (60) into (59), we obtain 


N 


r (T) = Aral Y\ 7 


1 


ti ko CO ] 4 

with the total phase $*0 (T) defined by 


sin 


A (T) 
* 


(61) 


$ko (T) = 


^fco (0) + Uko ( T ) 
2 



(62) 


Note that the total phase contains the large factor 1/7 ~ T, as expected. 
According to the estimates from the previous section, the minimal level sep- 
arations are estimated as Uko (0) ~ 1 and Uko (T) ~ Since in adiabatic 
limit the rate of parametric evolution is much smaller than the initial transi- 
tion frequency 7 << Uk 0 (0), the total phase is large and sin 2 CO] « 

Therefore, (61) reduces to 


N 


T{T) = 2 1 \ rlY. 


^ (Of 


(63) 


From (63), it follows that in order to evaluate the total failure probability 
of QAA, we need to know the density of states at the end of the evolution 
process. This is the consequence of the fact that the level separation is 
minimal at the end of the evolution. The final configuration of the optimal 
trajectory is given by (34) and corresponds to the standard GUE and the 
average density of states is given by the semicircle law (43). However, the 
expression (63) can not be evaluated in a continuous limit analogous to (38), 
since the corresponding integral diverges on the lower limit. This implies 
that the sum is accumulated on the lower limit corresponding to the low- 
frequency transitions, and the structure of DOS near the lower edge of the 
spectrum may be important. As it was shown in [14], the DOS in GUE near 
the ground state has the form of the power law 


p{\)~ \ q{x K, 


(64) 
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where the exponent q{ A) — 1/6 in the narrow range near the ground state 
and tends to the limit q (A) = 1/2 corresponding to (43), outside that range. 
Making use of (64), the probability (63) can be estimated as 


with 


N 


r (T) = 2-i 2 a 2 Q 


z! k 4 [Aw*] 


4 ’ 


(65) 


Ao> fc = (l/N) [piuj)}- 1 ~ (l/N) [Awfc] -9 k~ q . (66) 

Combining (65) and (66), we obtain 

r (T) = 2 7 V 0 2 N* f; 1 » ~t7 2(7 o ( 67 ) 

Jfe=l 

where p = 4/ (1 + 5) > 2. The estimate of the power of QAA is defined 
from the requirement that the failure probability (67) is small. Taking into 
account (53) and that in the GUE crfi ~ 1 /N, we estimate the runtime of the 
QAA 


T>>N ° • - 2 JTR)’ (68 > 

where the exponent q comes from (64). From (68), it follows that 17/14 < 
a < 5/6, where the upper and lower limits correspond to q — 1/6 and 
q = 1/2, respectively. The estimate (68) implies the exponential complexity 
of the QAA with a ~ 1. 

Comparing (46) and (67), we observe that the Markovian evolution model 
(12) may lead to polynomial complexity, whereas the driven RMT model 
(52) always leads to the exponential complexity of QAEA. This has a simple 
intuition. According to the above discussion, the main difference between 
the models (12) and (52) comes from different intertemporal level correla- 
tions. In the Markovian model (12), there are only short-range correlations, 
whereas in the driven RMT model (52), the long-range intertemporal corre- 
lations are present. On one hand, this leads to the complete cancellation of 
phases, which removes the exponentially small factors for the transitions to 
the highly excited states and leads to the diffusion type of dissipation process 
in (12). On the other hand, the baseness of the long-range intertemporal cor- 
relations leads to much weaker effective interaction between the states in the 
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Markovian model (12) as opposed to (52). The results (??) and (68) indicate 
that this ” weakness” of effective transitions is the leading effect ’that may in 
principle make the Markovian type QAEA (12) successful. 


5 Conclusion 

We analyze the power of quantum adiabatic evolution algorithms (QAEA) for 
solving random NP-hard optimization problems within a theoretical frame- 
work based on the random matrix theory (RMT). We present two types of 
the driven RMT models. 

In the first model, the ’’driving” Hamiltonian is represented by Brownian 
motion in the matrix space. We use the Brownian motion model t o obtain 
a description of multiple avoided crossing phenomena. The model allows for 
the close-form analytical treatment in adiabatic approximation, within the 
method of optimal trajectories [23]. We show that the failure mechanism of 
the QAA is due to the interaction of the ground state with the ’’cloud” formed 
by all the excited states. This confirms that in the driven RMT models, the 
Landau-Zener mechanism of dissipation based on the assumption of pairwise 
level interactions, is not important ([13]). We show that the QAEA may 
have a finite probability of success in a certain range of parameters, implying 
the possibility of polynomial complexity of the algorithm. This model can 
be viewed as a relatively ’’mild” test on the general QAA performance for 
the computationally hard optimization problems, in this case, the quantum 
evolution is not affected by the ” acciden tal” level crossing phenomenon that 
may lead to the QAA failure in standard models [11], and the only reason for 
the possible QAA failure is a complex structure of the problem Hamiltonian 
representing the cost function. 

The second model corresponds to the standard QAEA with the problem 
Hamiltonian taken from the Gaussian Unitary RMT ensemble (GUE). We 
show that the eigenvalues and eigenvector dynamics in this model can be 
mapped onto the dynamics in the Brownian motion model considered be- 
fore. This enables us to apply the optimal trajectories approach developed 
for the previous case, to the driven RMT model. However, the driven RMT 
and the Markovian Brownian m otion models have different structure of the 
intertemporal level correlations, and this leads to different performance of 
the QAEA. The Brownian motion model may lead to polynomial complex- 
ity, whereas the driven RMT model always gives the exponential complexity 



of the algorithm. In the Markovian Brownian motion model, there are only 
short-range correlations, whereas in the driven RMT model, the long-range 
intertemporal correlations are present. The absence of the long-range in- 
tertemporal correlations leads to much weaker effective interaction between 
the states in the Markovian model (12) as opposed to (52). Our results in- 
dicate that this ” weakness 7 ’ of effective transitions is the leading effect that 
may in principle make the Markovian type QAEA (12) successful. 



6 Appendix A 

According to (30) , the effective potential is given by 

U.tf ( A ) = - \ £ (vl (A) + ^7) , (69) 

with the velocities 1 % (A) defined by (22). From (22) and (69), we obtain 


^,(A) = 


k,l±k 


(70) 


( 7 & 

0* - * 1 ) 0* - V) A. (At - 


v ' '* ' v ' K ' ' l 7 fc,¥fc ( Afc ~ A '^ 

Under the cyclic permutation of the indices we have 


+ 


1 


+ 


1 


= 0 . 


(A* — A/) (A* — A i>) (A; - A//) (A; — A*) (A// — A*) (A v — A() 

therefore the first term in the square brackets in the r.h.s. of (70) reduces to 


E 


Wfitk (Afc - Ai) (A * ~ Ar) fcS ( Afc - A( ) 


(71) 


and the two terms in the square brackets in (70) sum up to zero. The second 
term in the r.h.s. of (70) is given by 


2 

1<T 


A, 2 A r (iV-l) 

— ' v/t 1 




— const . 


(72) 


Combining (70), (71) and (72) finally yields 


£ ? .//(A) = -^-E- x ‘- 


identical to (31) from the text. 


23 



7 Appendix B 

According to (32), the optimal trajectories satisfy 


A* (£) - i 2 \ k = 0. (73) 

The equations (32) are of the second order and therefore require two initial 
(or boundary) conditions for each k — 0 , ..JV. Suppose that we only have 
one initial condition for the optimal trajectories 

A* (t = 0) = A, (0) , (74) 

for each k and the final points A^ (t = T) are not fixed. Following the stan- 
dard procedure [9], w T e impose the additional ”transversality” boundary con- 
ditions at the final points 



which yield the boundary conditions for our problem 


(75) 


■ x. m - * [a mi - -7A, m + (g iprFMfj) ■ (76) 

The equations (32) with (74) and (76), provide a close set. From (32) with 
(74) and (76), we obtain 


A* (t) = A, (0) 


sinh [7 (T — £)] 


A fc (0 


sinh (7^) 


(77) 


sinh (7 T) sinh (7 T) 

The final ’’coordinates” {A*, (T)} are defined from (32) and (76) as solutions 
of the self-consistent equation 


A* (T) = X k (0) exp (-7 T) + 


[1 - exp (-27T)] 
7 



Afc (T) - X 3 (T) 


, ( 78 ) 


which indicates that in the large-time limit 7 T >> 1, the eigenvalues {A* (T)} 
are separated, i.e. the degeneracy is completely lifted, even though the initial 
distribution {A*, (0).} may have a significant degree of degeneracy. Specifically, 
(78) leads to the standard equilibrium GUE distribution described by (34) 
in the large-time limit. 
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8 Appendix C 

In this Appendix, we will derive the equilibrium iV-partiele distribution for 
the levels and show that it is indeed an equilibrium GUE distribution for 
any choice of the ’’rate” parameter a (t) introduced in the text. According 
to the first of equation (18), the non-stationary A- level distribution function 
P (A, t) satisfies the Fokker-Plank equation [18] 


^«^ + aX;[ v ‘ (A)j:>(A '‘ ) 

and the equilibrium distribution is given bv 


<7 2 dP (A, t) 
2 d\ k 


P eq (A) = C exp \-^V ef f (A) , (80) 

with a normalization constant C and the effective potential 

A 

V eq (A ) = ~ j dX k v k (A) = r X 2 k - a 1 ^ log |A fc - A,[ . (81) 

J k k,ljik 

Combining (80) and (81), we finally obtain 

P eq (A) = C exp lQ g (Afe - A;) 2 - 4 J2 A * » (82) 

k ^ k 

corresponding to the standard GUE distribution [7] for any choice of param- 
eters 7 and a. 


9 Appendix D 

In this Appendix, we employ the semicircle distribution to evaluate the prob- 
ability (38). According to (5), the levels of H (0) = Ho are equally spaced 
with u k o (0) = k and the degeneracy of the fc-th level is 

/n\ n\ , , 

9k = (jfc) = (n -*)!&!' ^ 

In the continuous limit of (43), we have 
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^ko (0) 


(84) 


= *(*/) = fa;, 

A o 

where the last equality follows from the observation that the width of the 
spectrum at t = 0 is 2u> n o (0) = 2 n and there are 2n+l distinct levels, whereas 
the spectrum at t — T is bounded within the range of 2A 0 . Combining (38) 
and (43), we obtain 


(P (T)) = — sinh (7T) 



(85) 


with the degeneracy gk( u ) given by (83) and (84). In the large n limit, the 
degeneracy (83) has a sharp peak at A; = nf 2, and (85) reduces to 


with 


(r(r)) = ^sinh (7 r)^ 



- w) 


Ap g jw) 

n lu ’ 


( 86 ) 


9(w) = V^ exp 

Evaluating (86) in the large n limit, we obtain 

(r (T)) = — — sinh (7T) (88) 

•7rn 7 Ao 

which reduces to (46) after taking (41) into account. 


n{u - A o) 2 

zxl 


10 Appendix E 

In this Appendix, we show that the estimate (??) for the power of the QAA 
algorithm given in the text for a particular case of the constant parameter 
characterizing the rate of changing the Hamiltonian a = const , can not be 
improved by introducing the time-dependent rate x (t)- 
The failure probability of the QAA is given by 
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2 

a 


( 89 ) 


N r t 

{T{t))=y / dt 1 s 

f^o (t')Y 


where u k Q (t) = \ k (t) — A 0 (t) axe the transition frequencies defined corre- 
sponding to the optimal trajectories {A fc (t)}. Under the transformation (??) 
with time-dependent evolution rate 7(f), the optimal trajectories satisfy, 
analogous to (32) 


A* (t) - 7 2 (t) A t = 0. (90) 

Analogous to (36), the solution of (90) can be presented in the form 

(t) ~ (0) yl (t) ~r Afc (T) (P 2 (t) , (91) 

where (<pi (f) , <p 2 (t)} are the two linearly-independent solutions of the second- 
order ODE (90) satisfying the following boundary conditions 


<pi( o) = l; <pi(T) = i; 

A2 (0) = 0; <p 2 (T) = l. 


(92) 


In this case, it is well known [24] that for the system (90), the Wronskian 


T/T/ ( ^ J 
r r y r j uicxi.jjLc;\a uy 


i-v (t) = (t) ^ (t) w 


dt 


dt 


= [^2 (f)i 2 4 


^1 (*) 


dt Lp 2 (t)\ ’ 


(93) 


is an integral of motion. Substituting (91) into (89), we obtain analogous to 
(37) 


(r CO) 


N .t 

= V / dfc- 


frj' -/o [w w (0) ¥>1 (t) + u > k 0 (T) <p 2 (*)f 

Making use of (93), (94) reduces to 


(94) 


<r (T)) 


du 


h'o w (y) [ll.%0 (0) U + Ulkd (U)] 


2 ’ 


(95) 
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where we made the substitution u = <pi ( t ) /</? 2 (f) and took (92) into account. 
After the integration (95) yields 


(T (T)> = 


^(7) 


E 


1 

u>ko (0) tdko (T) ’ 


(96) 


leading to the same bound on the runtime as the estimate (??). Note 
that for the case of constant x — 1/T considered in the text, the ba- 
sis functions {<pi (t) , y> 2 (£)} lead to (36), and the Wronskian is given by 
W (7) = 7/ sinh (7). In this case, (96) reduces to (38) from the text. 


11 Appendix F 

Following [17], consider the QAEA (2) 

H(t) = cos(f It) Ho + sin (Qt) V, (97) 

where Q is a real parameter characterizing the ’’speed” of the parametric 
evolution. The parametrization (97) is a particular case of (2) with 

a(t) = cos (Qt), ' (3 (t) = sin (fit) , T = ~r, (98) 

L 

and the parameters a(t) and (3 (t) being the monotonic functions of time 
during the runtime t e [0,T] of (97). From (51), we derive the parameters 
of the corresponding Ornstein-Uhlenbeck process in the form 


r (t) = ——In [cos(f2t)] , 
7 

- 27*o, 


. (99) 


implying that the trajectories (97) are mapped onto the Ornstein-Uhlenbeck 
process with constant parameters [17]. 
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